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1  Introduction 

Non-inva,sive  interrogating  techniques  are  most  valuable  in  determining  substructure  in  bi¬ 
ological  tissues  due  to  the  fact  that  they  usually  result  in  much  less  discomfort  in  subjects. 
Controlled  microwaves  (electromagnetic  waves  in  the  frequency  range  of  3  to  300  gigahertz) 
can  pass  through  many  media  without  causing  damage.  On  the  other  hand,  chemical  and 
physical  changes  in  biological  tissue  can  result  in  changes  in  its  electromagnetic  characteris¬ 
tics  such  as  electric  and  magnetic  polarization  mechanisms  and  conductivity.  Consequently, 
microwaves  sent  into  two  different  tissues  will  show  different  propagation  features,  and  anal¬ 
ysis  of  these  features  often  can  yield  useful  information  on  tissue  dysfunction.  Applications 
of  the  use  of  microwaves  in  non-invasive  interrogation  procedures  can  be  found  in  a  recently 
published  review  article  [1],  Use  of  ultrasonic  waves  is  another  popular  technique  used  in 
non-invasive  interrogation  of  media  in  both  industrial  and  medical  applications.  It  has  been 
well  known  since  1922  [5]  that  electromagnetic  and  sound  waves  can  interact  in  a  medium 
and  influence  each  other’s  propagation.  This  interaction  has  been  the  subject  of  substan¬ 
tial  investigation  in  acoustooptics  [7,  8,  10],  and  numerous  acoustooptic  devices  have  been 
developed  in  many  applications  in  industry  such  as  neural  nets,  optical  excision,  and  fiber 
optics  to  name  just  a  few. 

In  this  report,  we  investigate  the  possibility  of  using  the  interaction  between  electromag¬ 
netic  and  ultrasonic  waves  to  interrogate  the  structure  of  a  biological  medium.  The  medium 
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considered  here  is  a  fluid;  this  is  motivated  by  the  fact  that  the  major  component  of  human 
tissue  is  water.  We  focus  on  a  class  of  models  for  electric  polarization  in  the  context  of 
Maxwell’s  equations,  ft  is  well  known  (e.g.,  see  [2,  3]  and  the  references  therein),  that  po¬ 
larization  mechanisms  will  affect  substantially  the  propagation  of  an  electromagnetic  wave 
passing  through  a  medium.  Our  aim  here  is  to  demonstrate  how  modification  of  the  polar¬ 
ization  feature  in  a  tissue  by  an  ultrasonic  wave,  which  produces  a  virtual  interface,  can  be 
used.  An  electromagnetic  probe  sent  into  this  tissue  will  partially  reflect  from  this  artificial 
interface,  and  return  information  (e.g.,  geometry)  about  the  part  of  the  tissue  between  its 
surface  and  this  artificial  interface. 

For  the  propagation  of  the  microwaves,  we  assume,  as  in  [2,  3],  that  the  Maxwell’s 
equations  hold;  specifically, 


^  ^  <9B 

VxE  =  -7 5P 

(1.1) 

r)r> 

VXH  =  7-  +  J. 
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0) 
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V-B  =  0. 

(1.4) 

for  a  conductive  (Ohm’s  law)  dielectric: 

J  =  crE, 

(1.5) 

D  =  e0E  +  P, 

(1.6) 

B  =  goH  +  /i0M. 

(1.7) 

Here  E  is  the  electric  field  intensity,  D  is  the  (displacement)  electric  flux  density,  H  is  the 
magnetic  field  intensity,  B  is  the  magnetic  flux  intensity,  J  is  current  density,  pem  is  the  free 
(unpaired)  charge  density,  P  is  the  electric  polarization,  M  is  the  magnetic  polarization. 
In  our  initial  efforts,  we  concentrate  on  the  propagation  of  an  electromagnetic  wave  which 
is  uniform  in  the  x  —  y  plane  moving  in  the  2  direction.  The  motivation  (use  of  polarized 
impulsive  probes)  is  explained  more  fully  in  [2,  3].  This  assumption  allows  us  to  consider  the 
electromagnetic  fields  in  the  following  form:  E  =  E(f,z)i,P  =  P(t,z) i,  H  =  H(t,z) j.  We 
consider  electromagnetic  wave  propagation  in  the  normalized  interval  z  £  [0,  1],  assume  that 
the  fluid  slab  occupies  the  space  for  z  £  [zj,  1],  and  an  acoustic  wave  is  given  in  the  part  of 
the  fluid  for  z  £  [z2, 1]  with  0  <  Zi  <  z2  <  1.  Figure  1  is  a  sketch  of  the  geometry  considered 
here. 

Since  we  are  mainly  interested  in  biological  media  which  are  noil-magnetic,  we  further 
assume  that  the  magnetic  polarization  is  zero,  i.e.,  M  =  0.  Then  Maxwell’s  equations 
together  with  the  basic  constitutive  laws  yield  the  following  equation  in  the  domain  0  <  z  < 
1: 


d2E  a  dE  1  d2P  _  1  d2E  _  _  1  dJs 

dt 2  Co  dt  e0  dt 2  e0p0  dz 2  e0  dt 
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Figure  1:  Geometry  of  physical  problem. 


Applying  an  electronic  field  to  a  dielectric  material  will  cause  electric  polarization  which  is, 
of  course,  material  dependent.  In  our  investigation  here,  the  part  of  material  for  2  £  [0,  z-|], 
which  is  air,  is  assumed  to  have  zero  electric  polarization  and  zero  conductivity;  the  part 
of  material  for  2  £  [z1,z2]  is  assumed  to  obey  the  Debye  law  (Chapter  2,  [6]),  in  which  the 
electric  polarization  responds  to  the  electric  held  in  a  decaying  first  order  manner: 

dP 

t  +  P  =  eo(es  —  £oo  )E.  (1-9) 

Here  es  and  £<*,  are  the  static  relative  permittivity  and  high  frequency  relative  permittivity, 
respectively,  and  r  is  the  relaxation  time.  More  details  of  this  model  can  be  found  in  [6], 
and  other  models  for  the  polarization  can  be  found  in  [2,  3]  and  references  there. 

The  introduction  of  an  acoustic  wave  will  change  the  density  of  the  fluid  (indeed,  acoustic 
waves  are  simply  pressure  waves  which  involve  density  variations).  This  in  turn  will  affect 
electromagnetic  properties,  such  as  the  refraction  index,  of  the  fluid.  This  is  known  as  the 
acoustooptic  effect.  Consequently,  any  electronic  wave  transmitted  into  this  part  of  fluid  will 
be  modulated  by  the  acoustic  wave.  At  the  same  time,  the  material  electrostriction  caused 
by  the  electronic  waves  will  also  affect  the  propagation  of  the  pressure  wave  in  the  fluid  [9]. 
This  produces  a  fully  coupled  nonlinear  model  with  equations  for  both  the  electromagnetic 
and  acoustic  pressure  waves  (see  [4]  and  page  825  of  [9]).  In  our  initial  efforts,  we  focus  on 
the  effects  of  the  acoustic  wave  as  a  reflector  of  electromagnetic  waves.  We  ignore  the  effect 
of  electromagnetic  forces  in  the  acoustic  equation  under  the  tacit  assumption  that  the  effect 
is  weak.  To  demonstrate  the  effect  of  the  acoustic  wave  on  the  electronic  wave,  we  begin 
with  a  common  assumption  [8]  that  the  electric  susceptibility  is  an  afhne  function  of  the 


3 


acoustic  pressure  /;(/.,  z): 


X  =  Xo  +  Xi. P{t,z). 

Then  we  have 


P  =  =  £o(xo  +  XiPt1!  2))-E> 


(1.10) 


and  hence 


d2P  d2E  ,d2i> dpdE  d2E . 

r  -  eo^°-^F  +  e°^^7nT£  +  27i71i7+^li7T)- 


dt 


dt 2 


dt  dt 


dt2 


More  generally,  we  may  assume  that  the  fluid  in  the  acoustically  effected  part  of  the  domain 
obeys  a  generalized  pressure  dependent  polarization  rule  (Chapter  9,  [6]): 

1  d2P  x/  dE  £.  d2E 

7„W  =  Mp)e  +  Mp)TT  +  Mp)W 


To  simplify  issues  in  this  preliminary  investigation,  we  take 


fo{p)  =  0,  fi(p)  —  0,  f2(p)  =  Kp(t,z), 


which  we  note  is  not  a  special  case  of  (1.10).  Thus  our  polarization  assumption  in  the 
acoustic  fluid  is  given  by 


1  d2P  .  d2E 

70lF  =  Kpltz)lF' 


(in) 


Using  (1.11)  in  (1.8),  we  have  the  following  partial  differential  equation  for  the  electric  held 
in  the  region  disturbed  by  acoustic  waves: 


d2E  a  dE  .  d2E  1  d2E  1  dJs 

-aF  +  77r  +  Kp{t’z)-aF  =  7i777-77r-  Z2<Z<1- 

This  equation  is  very  similar  to  the  one  given  in  [9]  derived  from  thermodynamical  consid¬ 
erations. 

To  complete  the  demonstration  model,  we  assume  that  the  material  outside  the  space 
%  (E  [0, 1]  can  absorb  the  electronic  wave  completely;  hence  we  can  use  the  following  boundary 
conditions  for  the  electric  held: 


with 


dE,  „  dE,  , 
^(b°)  -c— (f,°)  -  0, 

dE,  ^  c  dE 

-rr-(b  1)  H - / 

yjl  +  Kp[t:  Z  ) 


(M)  =  0, 


Co/^'O 
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We  also  assume  the  initial  conditions: 


£(0,z)  =  0,  —  (0,z)  =  0. 

For  the  purpose  of  deriving  an  efficient  numerical  scheme,  we  reduce  the  second  order 
derivative  of  P  in  (1.8)  by  the  Debye’s  polarization  law  (1.9).  This  results  in  a  first  order 
derivative  of  P  along  with  a  similar  additional  term  for  E  in  the  basic  Maxwell’s  equation. 
Thus,  we  use  the  following  initial-boundary  value  problem  to  model  the  dynamics  of  the 
electromagnetic  fields: 


d2E  dE  dP 

“('■*)«>+*(*)« +'W  a 

d2E 

=  d(z)  +  E(t,z),  z  G  (0, 1), 

(1.12) 

dP 

-jjj-  —  72 P  +  I2E,  z  €  [zi,  z 2], 

(1.13) 

dE,  ,  /—~dE1 

dt  M)  yfd/a  {t,  0)  =  0, 

dE,  ,  rr-dE , 

"^"(t  1)  +  yd/a—(t ,  1)  -  0, 

(1.14) 

O 
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O 
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II 
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(1.15) 
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(1.16) 

where 


a(t,z)  =  1  +  k  p(t,z)x[Z2,i], 

H'-)  —  Xffi,1]  “I 

£o  eo 

e\Z)  X[zi,z2]l 
Co 

d(z)  =  c2  =  - , 

eo/^o 


1  ^o(Cs  ^OO  )  771/1  \  1  9Js 

m  = - ,  72  =  - F(t,  z)  = - 

t  t  e0  at 


and  X[a?i ,ora]  is  the  usual  characteristic  function  for  the  interval  [.7q ,  .r^J . 


2  A  finite  element  scheme 

To  approximate  the  solutions  to  the  unknowns  /•.  (/.,  z)  and  P(t,  z)  in  the  model  presented  in 
the  previous  section,  we  hrst  introduce  a  partition  in  the  space  variable: 

0  =  Xq  <  x1  <  ■  ■  ■  <  xNz  =  1, 


such  that 
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for  some  integers  1  <  i\  <  i2  <  Nz,  and 


h  =  max  (x,-  —  x.-i ). 

Then  we  let  Sh  C  H1(  0, 1)  be  the  standard  linear  finite  element  space  defined  on  this  partition 
with  a  set  of  basis  functions  We  also  use  a  partition  in  the  time  variable: 

0  =  t°  <  t1  <  ■  ■  ■  <  tNt  =  T, 

with  A t  =  tl  —  f*_1,  i  =  1,  •  •  ■  ,  Nt. 

As  usual,  we  start  with  the  weak  form  (see  [2,  3])  of  the  model  partial  differential  equation 

(1.12): 

d^E  dE 

(vi  a~^r)  +  ( v >  h~fr)  +  (a  er/2  p)  +  (a  e72  e) 

=  I  hi  _  r  -  +  A  A 

for  any  v  (E  ih1(0,l),  where  (  ,  )  denotes  the  usual  T2(0,  1)  inner  product.  From  this,  we 
can  formulate  the  following  finite  element  discretization  of  (1.12): 

(4  ad»En)  +  (&,  IT)  +  (<4  er/2  Pn)  +  (<4  e72£n’1/4)  =  -V^ddtEnfo\z=1 

-\^ddtEn<j>t\z=°  -  (^,4 E*#4)')  +  (^l#F(fV)),  (2.17) 

where 


N 

i=i 

N 

Pn(z)  =  '£p]Mz)*P(tn,z), 

3  =  1 


and  the  following  finite  difference  notation  is  used  for  the  average 

[fn-H  +  2ff,  +  cm-1 

4 

along  with  the  difference  quotients 


dttUn  = 


f7”+1  -  2f7”  +  [A”1 
A?  ' 


at/n  = 


[/n+l  _  [fn-l 

2A/ 


As  a  scheme  to  compute  /A,  assuming  that  Pn  is  given  exactly,  we  expect  the  above  scheme 
to  have  the  following  accuracy; 


\\En  -  E(tn,-)\\L2  <  C{h2  +  At2), 
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provided  that  the  solution  is  smooth  enough.  To  update  lr"  at  each  time  step,  it  is  preferable 
to  choose  a  scheme  that  has  good  stability  and  whose  accuracy  matches  that  used  to  compute 
En .  Hence  we  use  the  following  A- stable  second  order  Adams-Moulton  scheme  to  discretize 
the  equation  (1.13): 


pn 


_Lm  (p»-i  +  A(72(£*  +  e-')+  . 


(2.18) 


According  to  the  given  initial  conditions,  we  should  set 


E°  =  0,  P°  =  0, 


(2.19) 


for  the  finite  element  approximations  at  the  0-tli  time  level.  By  the  Taylor  expansion  and 
applying  (1.12),  (1.13),  (1.15),  and  (1.16),  we  have 


dt 2 


At2 


2a(0,  z) 


F(  0,4 


Therefore  we  can  let 


E1  = 


At 2 


2a ( 0, z) 


F(  0, 


(2.20) 


Putting  all  of  these  discretizations  together,  we  have  the  following  algorithm  to  generate 
approximations  to  both  E{t,z)  and  l>(l..z): 

Step  1  .  Compute  E°  and  P°  by  (2.19). 

Step  2  .  Compute  E1  by  (2.20). 

Step  3  .  Then  for  each  n  =  1, 2,  ■  ■  ■  ,  Nt  —  1,  we  use  (2.18)  to  compute  Pn  ~  P(tn,  z)  and  use 
(2.17)  to  compute  En+1  fa  E{tn+1,z). 


3  Some  numerical  simulations 

To  test  the  model,  we  assume  that  a  time  “windowed”  electromagnetic  point  source  input 
(e.g.,  see  [2,  3])  is  given  at  the  left  boundary  point  z  —  0  such  that 

■Js{t,z)  =  -S(z)x[o,tf]{t)sin(tds  /.). 

The  frequency  in  the  source  is  assumed  to  be  in  the  microwave  range,  i.e. , 

G  [3  x  109 Hz, 3  x  1011  Hz]. 
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The  pressure  is  given  by 


K,p(l,.  z)  =  ng{z)sin\LOvt), 

such  that  its  frequency  is  in  the  ultrasonic  range,  i.e. , 

u>p  £  [0.lMHz,2hMHz\, 

and  in  all  the  computations  presented  here,  g(z)  =  X[z2,i](z)' 
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^oo 
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Up 

7T  X  10  Y 

Us 

7T  X  101U 

At 

tf/ 1600 

h 

1/900 

Table  1:  Parameters  used  in  Example  1. 

Example  1:  We  chose  all  the  parameters  in  this  example  similar  to  those  for  water  except 
for  the  relaxation  time  r,  with  the  value  for  r  used  here  somewhat  larger  than  those  given  in 
the  literature.  We  observed  that  smaller  r  leads  to  a  weaker  transmitted  electronic  wave  into 
the  fluid  which  hinders  observation  of  the  interaction  between  the  electronic  and  acoustic 
waves.  Some  of  the  parameters  are  listed  in  Table  1.  Representative  plots  are  in  Figure  2-5. 

Example  2:  The  data  used  in  this  example  are  listed  in  Table  2.  All  the  values  are  the  same 
as  those  in  the  previous  example  except  for  the  frequency  in  the  acoustic  wave.  In  this  case, 
a  lower  frequency  is  used  in  the  acoustic  wave  and  we  notice  that  the  reflected  wave  from 
the  acoustic  beam  is  weaker  than  that  in  the  previous  example,  see  Figure  6. 


tf 

4.0x  10_1U 
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^OO 
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u 
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Table  2:  Parameters  used  for  Example  2. 
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time  level  11001 


Figure  4:  Example  1.  The  electronic  wave  in  the  part  of  fluid  [z2, 1]  where  there  is  an  acoustic 
wave.  A  reflected  wave  generated  from  the  acoustic  beam  is  seen  in  [zi,z2]- 


time  level  17001 


Figure  5:  Example  1.  The  reflected  wave  from  the  acoustic  beam  approaches  the  left  bound¬ 
ary. 
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time  level  11001 


Figure  6:  Example  2.  The  electronic  wave  in  the  part  of  the  fluid  where  there  is  an  acoustic 
wave.  A  rather  weak  reflected  wave  is  generated  from  the  acoustic  beam  at  a  lower  frequency. 


Example  3:  A  smaller  relaxation  time  r  is  used  in  this  example,  all  other  parameters  are 
the  same  those  in  Example  1  (see  Table  3).  From  Figure  7,  8,  and  9  we  can  see  that  the 
transmitted  electronic  wave  into  the  fluid  is  weak,  with  most  of  the  electronic  wave  reflected 
from  the  interface  between  air  and  fluid. 


tf 

4.0x  10_1U 
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8.1  x  10-12 

eoo 

5.5 

£ s 

78.2 

u 

l.oxicr5 

K 
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7 r  x  10v 

Us 
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At 

tfl 1600 
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1/900 

Table  3:  Parameters  used  for  example  3. 
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time  level  6551 


Figure  9:  Example  3.  Most  of  the  wave  is  reflected  back  into  the  domain  filled  with  air. 

4  Identification  of  material  depth 

In  this  section,  we  consider  the  problem  of  using  the  signal  E(t,  0)  collected  at  the  left 
boundary  z  =  0  to  estimate  the  material  depth.  To  be  specific,  we  assume  that  all  the 
physical  parameters  other  than  the  depth  of  the  material  are  given  (these  can  be  estimated 
using  the  first  reflected  waves  from  the  interface  at  z1,  see  [2,  3]),  and  the  interface  of  the 
material  closer  to  the  source  (which  is  also  at  z  =  0)  is  fixed  at  'z  =  Z\.  As  before,  we  denote 
the  position  of  the  interior  interface  between  the  fluid  and  acoustic  domain  by  z2.  After  an 
electromagnetic  wave  is  sent  to  the  material,  reflected  waves  will  be  generated  at  the  two 
interfaces  of  the  material  that  will  propagate  back  to  the  left  boundary  where  the  wave  is 
generated.  Intuitively,  the  difference  between  the  times  when  the  wave  reflected  from  the 
first  interface  and  that  from  the  second  interface  reach  the  left  boundary  depends  on  the 
depth.  Hence  we  shall  attempt  to  estimate  the  depth  of  the  material  from  this  difference. 

We  first  need  a  procedure  to  detect  the  time  when  a  wave  reflected  from  the  interface 
reaches  the  left  boundary.  Figure  10  depicts  a  typical  data  function  /•:  (/.,  0).  Note  that 
the  data  function  is  essentially  zero  except  in  three  subintervals.  We  first  considered  using 
the  velocity  Et(t,  0)  to  estimate  z2.  We  expect  the  function  Et(t,  0)  to  become  zero  after 
the  source  has  been  turned  off  for  a  while.  Then  Et(t,  0)  becomes  nonzero  when  the  wave 
reflected  from  the  first  interface  arrives  at  the  left  boundary.  This  is  followed  by  a  period 
of  time  when  Et(t,  0)  becomes  zero  again.  Afterwards,  Et(t,  0)  becomes  nonzero  due  to  the 
arrival  of  the  wave  reflected  from  the  second  interface.  Figure  11  presents  a  typical  plot  of 
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0)  |.  To  aid  in  our  estimation  procedure,  we  define  Ty  to  be  the  first  time  when 

/v,(7i,0)  >  Ty  >  C,.. 

where  both  constants  Cy  and  C2  can  be  determined  from  the  measurement  of  E(t,  0).  The 
time  C2  should  be  larger  than  the  time  within  which  the  source  is  generated.  The  constant  Cy 
should  be  chosen  large  enough  to  distinguish  the  change  in  data  due  to  measurement  error. 
Since  Zy  is  assumed  to  be  fixed  in  the  identification  of  depth,  Ty  should  be  independent  of 
z2.  We  next  let  T2  be  the  first  time  such  that 

/'  ;•(  /2,  0)  >  C3,  T2  >  Ty  +  C4. 

The  constant  C3  can  be  chosen  in  a  way  similar  to  that  for  Cy ,  which  we  can  simply  take 
C3  =  Cy.  The  value  of  64  should  be  chosen  so  that  the  wave  reflected  from  the  first  interface 
has  passed  the  left  boundary  at  the  time  /  =  /2  +  C4.  Clearly  T2  depends  on  the  position 
of  the  second  interface,  and  we  denote  it  by  T2  =  7 2 ( z2 )  .  Then  we  define  a  function  of  z2  as 
the  difference  of  these  two  characteristic  times: 

L(Z2)  =  T2(Z2)-Ty, 

and  our  identification  problem  leads  to  looking  for  a  value  z2  such  that 

L(z2)  =  L(z-)t  (4.21) 

where  z2  is  the  true  position  of  the  second  interface  of  the  material  from  which  the  data,  is 
collected.  Note  that  ^2  is  unknown,  but  L(z2)  can  be  generated  from  the  data. 

To  see  the  behavior  of  T(z2),  we  calculated  its  values  for  various  z2  in  the  neighborhood 
of  z2.  A  typical  plot  of  T(^2),  given  in  Figure  12,  suggests  that  L(z2)  acts  almost  like  an 
affine  function  in  the  neighborhood  of  z2.  In  this  case,  the  secant  method  is  a  good  candidate 
for  computing  z2  from  (4.21). 

Example  4:  We  assume  that  the  true  position  of  the  second  interface  is  at  z2  =  2/3.  The 
data  E(t,  0)  was  generated  by  the  finite  element  scheme  presented  in  Section  2  with  z2  =  z*2 
and  other  parameters  listed  in  Table  4.  The  plot  for  the  data  of  E(t ,  0)  is  in  Figure  10,  and 
the  absolute  value  of  its  numerical  derivative  is  plotted  in  Figure  11.  From  these  plots,  we 
decided  to  choose  the  constants  Cy,  C2,  C3:  C4  as  follows: 

Cy  =  C3  =  4.5  x  1010,  C2  =  2.0  x  10“9,  C4  =  1.25  x  10“9. 

With  these  constants,  we  can  find  that 

L(z*2 )  =  2.2165  x  10“9. 

Then  we  used  the  secant  method  to  solve  iteratively  for  z2  in 

L(z2)  =  2.2165  x  10“9, 
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2.35 


Figure  12:  A  typical  plot  of  Llyz2)  in  the  neighborhood  of  z2. 
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Table  4:  Parameters  used  for  Example  4. 


and  obtained  the  following  approximation  to  the  exact  position  of  the  second  interface: 

i2  =  0.66653017606683  »  z2*  =  2/3. 

We  remark  that  in  the  computation  of  z2,  we  used  h  =  1/900  and  r  =  tf/ 1600  in  the  initial 
boundary  value  problem. 

We  observe  that  the  above  scheme  appears  to  be  a  reasonable  approach  only  if  we  can  have 
a  dependable  measurement  for  the  velocity  Et(t,  0).  This  approach  most  likely  will  not  work 
if  numerical  differentiation  must  be  used  to  obtain  an  estimate  of  Et(t,  0)  from  a  measurement 
of  E(t ,  0)  with  noise.  This  is  due  to  the  well  known  catastrophic  behavior  encountered  in 
using  numerical  differentiation  on  error-polluted  data.  Even  a  very  small  amount  of  error  in 
the  data  for  E(t ,  0)  will  make  the  estimate  of  /'/<(/..  0)  generated  by  numerical  differentiation 
meaningless.  For  example,  if  we  pollute  the  data  given  in  the  previous  example  by  a  uniformly 
distributed  random  relative  error  with  a  magnitude  only  5%  of  that  of  the  data,  then  the 
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Figure  13:  A  typical  plot  of  an  estimated  A  f  i.  0  )  |  generated  from  numerical  differentiation 
on  a  measurement  of  K( I  .  0)  with  a  random  error  whose  magnitude  is  only  5%  of  the  true 
data. 

Et(t,  0)  plotted  in  Figure  13  generated  by  numerical  differentiation  does  not  give  any  useful 
information  about  the  time  lag  between  the  two  reflected  waves. 

On  the  other  hand,  the  pulse  signal  E(t,  0)  itself  appears  to  be  rather  robust  with  respect 
to  the  random  noise  from  the  point  of  view  of  indicating  the  time  lag  between  the  reflected 
waves  from  the  two  interfaces.  Figure  14  is  a  plot  of  the  absolute  value  of  the  data  for 
E(t,  0).  Adding  a  uniformly  distributed  random  error  at  the  5%  level  yields  data  plotted 
in  Figure  15  from  which  we  can  still  easily  discern  the  times  when  the  two  reflected  waves' 
arrive  at  the  left  boundary.  Hence  we  introduce  another  quantity  to  describe  the  time  lag 
from  a  measurement  of  E(t ,  0)  as  follows. 

We  define  'l\  to  be  the  first  time  such  that 

£(Ti,0)  >  Cu  Tj>C2 , 

where  both  constants  C\  and  C2  can  be  determined  from  the  measurement  of  E(t,  0).  The 
value  of  C2  should  be  larger  than  the  time  during  which  the  source  is  generated.  The 
constant  C\  should  be  chosen  sufficiently  large  so  as  to  distinguish  the  change  in  data  due 
to  measurement  error.  The  time  7\  should  be  independent  of  z2  since  Zj  is  assumed  to  be 
fixed  in  the  estimation  of  depth.  We  let  T2  be  the  first  time  such  that 

£(T2,0)>A3,  72  >  7i  +  C4. 
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The  constant  C3  can  be  chosen  in  a  way  similar  to  that  for  C\  (we  can  even  simply  let 
C3  =  Ci).  The  constant  C4  should  be  chosen  so  that  the  wave  reflected  from  the  first 
interface  has  passed  the  left  boundary  at  the  time  t  =  T2  +  C4.  As  expected  T2  depends 
011  the  position  of  the  second  interface,  and  we  denote  it  by  T2  =  T2{z2).  Then  we  define  a 
function  of  z2  as  the  difference  of  these  two  characteristic  times: 

L{z2)  =  T2(z2)-T1. 

Obviously,  the  construction  procedure  for  L(z2 )  is  similar  to  the  one  associated  with  data 
Et{t ,  0),  but  no  differentiation  is  used.  As  before,  we  define  the  solution  to  the  identification 
problem  as  a  quantity  z2  that  satisfies 

L(z2)  =  L(z*2).  (4.22) 

Since  this  new  time  lag  function  of  z2  also  appears  to  behave  linearly  in  the  neighborhood 
of  the  exact  location  of  the  second  interface  z2  (see  Figure  16),  we  again  believe  that  the 
secant  method  is  a  good  candidate  for  computing  z2  from  (4.22). 

Example  5:  The  function  L[z2)  just  introduced  is  also  rather  insensitive  to  noise  in  the  data. 
See  Table  5  for  its  behavior  with  respect  to  the  random  error  with  a  uniform  distribution, 
and  Table  6  for  its  behavior  with  respect  to  the  random  error  with  a  normal  distribution. 


Noise  Level 

L(z2)  by  datum  with  noise 

0 

2.223  x  10“9 

1% 

2.22325  x  10“9 

5% 

2.2225  x  10“9 

10% 

2.2225  x  10“9 

15% 

2.22275  x  10“9 

20% 

2.22225  x  10“9 

Table  5:  Measure  of  sensitivity  of  L[z2)  to  the  noise  with  a  uniform  distribution  in  the  data. 

Here  we  used  the  finite  element  solution  to  generate  a  data  for  E{t ,  0)  with  z*2  =  2/3  and 
other  parameters  listed  in  Table  4.  The  plot  for  the  data  of  E(t,  0)  is  in  Figure  10,  and  the 
absolute  value  of  the  two  reflected  waves  received  at  the  left  boundary  is  plotted  in  Figure 
14.  From  these  plots,  we  decide  to  choose  the  constants  Cj,  C2,  C3,  C4  for  the  definition  of 
L[z2)  as  follows: 

Ci  =  C3  =  1.0,  C2  =  2.0  x  10“9,  C4  =  1.75“9. 

We  perturbed  this  data  by  random  numbers  with  various  noise  level,  and  used  the  error 
polluted  data  to  generate  values  of  L[z2)  in  Table  5  and  Table  6.  Note  that  L(z2)  changes 
little  even  with  the  data  at  the  20%  noise  level. 
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Standard  deviation 

L[z2)  by  datum  with  noise 

0 

2.223  x  10“y 

0.01 

2.2235  x  10“9 

0.05 

2.222  x  10“9 

0.1 

2.2203  x  10“y 

0.15 

2.2258  x  10“y 

0.2 

2.219  x  10“y 

Table  6:  Measure  of  sensitivity  of  L[z2)  to  the  noise  with  a  normal  distribution  in  the  data. 


Example  6:  The  estimation  procedure  for  the  depth  z*2  based  on  this  new  time  lag  function 
works  well  with  data  polluted  with  random  errors  at  various  levels  as  seen  Table  7.  Even 
the  data,  with  10%  noise  yields  a  good  approximation  to  the  location  of  the  second  interface. 
Note  that  in  the  computation  of  z2.  we  used  the  physical  parameters  listed  in  Table  4,  but 
we  used  h  =  1/900  and  At  =  tfj  1600  to  solve  the  initial  boundary  value  problem,  and  the 
exact  location  of  the  second  interface  is  z*2  —  2/3. 


Noise  Level 

z2  by  data  with  noise 

0 

0.66632745773773 

1% 

0.66632745773773 

5% 

0.66640117743327 

10% 

0.66629059787037 

Table  7:  The  estimation  procedure  is  rather  robust  with  respect  to  the  random  noise  in  the 
data.  Uniformly  distributed  random  errors  were  used  in  these  computations. 
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2.35 


Figure  16:  A  typical  plot  of  i(z2)  in  the  neighborhood  of  z\  —  2/3. 

5  Concluding  remarks 

In  the  preliminary  investigations  reported  on  in  this  note,  we  have  demonstrated  the  potential 
to  employ  internal  acoustic  fields  as  reflectors  for  electromagnetic  probes  in  the  interrogation 
of  dielectric  media.  The  associated  inverse  problems  are  based  on  time  domain  formulation 
of  the  acoustooptic  signals  in  the  media.  Encouraged  by  these  early  findings,  our  efforts 
are  continuing  with  more  involved  electric  polarization  models  as  well  as  with  models  of  the 
acoustooptic  interaction  that  are  closer  to  physical  reality. 
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